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We consider dark matter consisting of weakly interacting massive particles (WIMPs) and revisit 
in detail its thermal evolution in the early universe, with a particular focus on models where the 
annihilation rate is enhanced by the Sommerfeld effect. After chemical decoupling, or freeze-out, 
dark matter no longer annihilates but is still kept in local thermal equilibrium due to scattering 
events with the much more abundant standard model particles. During kinetic decoupling, even 
these processes stop to be effective, which eventually sets the scale for a small-scale cutoff in the 
matter density fluctuations. Afterwards, the WIMP temperature decreases more quickly than the 
heat bath temperature, which causes dark matter to reenter an era of annihilation if the cross-section 
is enhanced by the Sommerfeld effect. Ifere, we give a detailed and self-consistent description of 
these effects. As an application, we consider the phenomenology of simple leptophilic models that 
have been discussed in the literature and find that the relic abundance can be affected by as much 
two orders of magnitude or more. We also compute the mass of the smallest dark matter subhalos 
in these models and find it to be in the range 0{1G-^°Mq) - 0(10 M©); even much larger cutoff 
values are possible if the WIMPs couple to force carriers lighter than about 100 MeV. We point out 
that a precise determination of the cutoff mass allows to infer new limits on the model parameters, 
in particular from gamma-ray observations of galaxy clusters, that are highly complementary to 
existing constraints from g — 2 or beam dump experiments. 

PACS numbers: 95.35. +d, 98.80.-k, 98.80.Cq, 95.30.Tg 



I. INTRODUCTION 

A plethora of independent cosmological observations, 
taken over a large range of distance scales, strongly sup- 
ports the existence of cold, non-baryonic dark matter 
(DM); according to the most recent estimates, it con- 
tributes a fraction of SIdm = 0.229 ± 0.015 to the total 
energy content of the universe [T]. Weakly interacting 
massive particles (WIMPs) are the leading candidate for 
the so far still obscure nature of the DM, well motivated 
from particle physics and potentially offering promising 
prospects for detection in indirect, direct and accelerator 
searches [2]. When mentioning DM in this article, we 
will always refer to WIMPs. 

One of the most appealing aspects of WIMPs is their 
thermal production in the early universe, resulting in a 
relic density that is roughly of the same order as the ob- 
served DM density today [3 . While this somewhat sim- 
plistic picture actually changes when taking into account 
complications ^ like co-annihilations, thresholds or res- 
onances, the main idea still remains the same: at very 
high temperatures, WIMPs are kept in thermal equilib- 
rium through annihilation and creation processes with 
the standard model particles of the heat bath; once the 
interaction rate falls considerably behind the Hubble ex- 
pansion, the WIMP number density 'freezes out' and is 
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only affected by the further expansion of the universe (for 
a very thorough treatment, see e.g. Refs. [5, 6 ). 

Even after this chemical decoupling, WIMPs stay in 
thermal contact with the heat bath due to the much 
more frequent elastic scattering events with the relativis- 
tic standard model particles [3 H]. At kinetic decou- 
pling, which actually happens on a rather short timescale 
[9l [To] , even these scattering processes cease to be effec- 
tive; in the standard picture, there is then no further in- 
teraction between WIMPs and the heat bath. From this 
point on, density perturbations in the DM component 
are thus no longer tightly bound to density perturbations 
in the radiation component and can freely develop, with 
free streaming effects setting a lower bound to the size 
of over-dense regions that will eventually collapse under 
the influence of gravitation [TTHT3] ; acoustic oscillations 
lead to a similar cutoff [121 [TS] . 

While the processes that determine the relic density 
and the cutoff in the power spectrum of matter den- 
sity fluctuations often happen at completely unrelated 
timescales, however, this does not necessarily have to be 
the case. In models where the annihilation rate at small 
relative velocities is strongly enhanced, e.g., the effects of 
annihilation and scattering processes on the DM distri- 
bution can strongly interfere and change the simple pic- 
ture sketched above. In particular, one may encounter 
a new era of DM annihilations after kinetic decoupling 
has taken place [m - fTOj . Physically, such strongly en- 
hanced annihilation rates can be motivated by the re- 
peated exchange of relatively light force carriers, known 
as the Sommerfeld effect [20] ; in the context of DM anni- 
hilation, this has first been studied for heavy neutralino 
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DM coupling to SU (2) gauge bosons pTH25] , but seen a 
boost of interest (see, e.g. Refs. [24H30] ) in the context 
of possible DM explanations for the apparent anomalies 
in the observed cosmic ray electron and positron fluxes 

Here, we revisit in detail how WIMPs decouple from 
the thermal bath in the early universe and derive, by 
extending the standard calculation, a coupled system of 
equations that describes the evolution of the WIMP num- 
ber density and velocity dispersion (aka the WIMP "tem- 
perature" ) . This allows to compute the evolution of these 
quantities to a very high precision and, in particular, 
treat situations like the one mentioned in the previous 
paragraph in a fully consistent way that can be applied 
to any DM model where the annihilation rate is enhanced 
for small relative velocities. 

As an application, we focus on a simple toy-model with 
large Sommerfeld enhancement and compute both the 
final relic density and the size of the smallest subhalo 
masses. We demonstrate that the above-mentioned ef- 
fects can be sizable and that the consistent treatment de- 
rived here is indeed necessary to reliably compute these 
quantities in DM models where the Sommerfeld effect 
is relevant. Interestingly, gamma rays from DM annihi- 
lation in galaxy clusters may place strong lower limits 
on the cutoff mass in this kind of models |32l |33]; we 
demonstrate here that a reliable and self-consistent com- 
putation of this quantity can therefore be used, in prin- 
ciple, to translate these limits into constraints on the pa- 
rameter space that are highly complementary to existing 
constraints from g — 2 or beam dump experiments. 

This article is organized as follows. We start in Section 
ITTl with a concise review of chemical and kinetic decou- 
pling of WIMPs, extending the standard treatment in 
such a way as to take into account possible interferences 
between the two. In Section III wc then discuss in detail 



the evolution of the DM phase-space distribution after ki- 
netic decoupling, with a particular focus on the possibil- 
ity of Sommerfeld-enhanced annihilation rates. Section 
IV is devoted to the application of our general formalism 
to a concrete leptophilic DM particle model. We discuss 
our results in Section [V| and present our conclusions and 
an outlook in Section fVll In Appendix [X} we collect some 
relevant details about the Sommerfeld effect. For the lep- 
tophilic toy model studied in Section |TVj we provide cur- 
rent constraints in Appendix [B] and, for convenience, a 
collection of scattering and annihilation matrix elements 
in Appendix [C| 



II. THE STANDARD THERMAL EVOLUTION 
OF WIMPS 

In this section, we review the standard case of ther- 
mally produced particles in the expanding universe. The 
evolution of their phase-space density /(p), in partic- 
ular, is described by the Boltzmann equation which in 
a Friedmann- Robertson- Walker metric reads (see, e.g.. 



E [dt - i/p • Vp) / = C[f] . 



(1) 



Here, p'' = {E, p) denote comoving WIMP momenta 
and H — a /a the Hubble parameter; all interactions 
are contained in the collision term on the right-hand 
side. As we will see in detail, this equation describes 
how WIMPs start in perfect thermodynamic equilibrium 
with the early, very dense and hot universe and eventu- 
ally decouple completely from the heat bath, in several 
distinct stages, as the universe continues to expand. 



A. Chemical decoupling 

The evolution of the DM particle number density 



9x 



;/(p) 



(2) 



is affected both by the expansion of the universe and by 
DM annihilation into, or creation from, SM particles. Re- 
stricting ourselves to 2-body processes, the corresponding 
contribution to the collision term is given by 



25x 



^ , (27r)32a;7 {2T:f2Cj J (2Ti f2E 
x{2-Kf5^^\p + p-k-k) 



9xEJ2 J (^"-i'^: 

X [/eq(£;)/cq(^) - f{E)f{E) 



(3) 



where fc^ = (w,k) and fc^ — (w,k) are the 4-momenta of 
the SM particles X and g = ffcq = (e'^/'^ tti'^ir 
distribution functions; |A^|^ refers to the matrix ele- 
ment summed over both SM and DM internal (spin) de- 
grees of freedom gx.x- Note that we can safely neglect 
Pauli-suppression factors in the non-relativistic regime 
((P^) ^ "7,^) we are considering here. The second step 
follows from CP invariance and the fact that annihilation 
and creation processes should happen with the same fre- 
quency in equilibrium; the velocity appearing here is the 

M0ller velocity, v^d = vm&i = {EEy^J{p ■ pY - m\. 



In order to proceed, one usually assumes that f{E) cx 
fcqiE) = e~^/"^, with a factor of proportionality that de- 
scribes an effective chemical potential which may depend 
on T (but not E); this is motivated by the fact that the 
much more abundant scattering processes of DM with 
SM particles still keep the DM particles in kinetic, but 
not chemical equilibrium - see also the following Section. 
Integrating Eq. ([ij over J cPpg^/ [(27r)'^£'] then results 
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+ 3Hn^ = -(crWroOcq [nl - n°^'^ 



(4) 



with n'^ = ra^g^K2{x) / {2'r:'^x) in the non-relativistic 
regime and 

/ \ = A_ f Jlr. f JlL 

^'^^-'^'^'^ - „cq2 J (2^)3 J (2^)3 

XVrclCTxx^-^A-/cq(£^)/cq(^) (5) 

4x3/2 n 



(6) 



where K2 is the modified Bessel function of second order 
and V is the velocity of each WIMP in the center-of-mass 
frame. The last approximation is valid for large x and, 
for X > 10, reproduces the full analytical result to 
an accuracy of better than 1% for all relevant functional 
dependences of (av^Qi). A simple power-law 



in particular, results in 

, , 2(To ^ / 3\ _ _ _ 

(o-l^rcOcq ^ r\^n + -j X = UqX 



(7) 



(8) 



The next step is to introduce dimensionless variables 



X = m^/T , 
Y = n^/s, 



(9) 
(10) 



where T is the temperature of the heat bath and s 
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g*s{T)'^T^ the entropy density. Using entropy conser- 



vation, dt (a^s) = 0, one can then transform Eq. Q to 
the convenient form 



Y 



V Sg.sJ Hx I r2 ) ' ^^'^ 



where ' = d/dx. Note that without the assumption 
of f{E) (X /eq(i?), we would have arrived at a very 
similar equation, with (l — Y^^/Y^) being replaced by 

((cr^;rei)/(o-?;rei)eq_- Y^JY'^), where (avrei) is defined ex- 
actly as in Eq. pi), but for an arbitrary WIMP distri- 
bution, i.e. with (/eq, 71^*5 ) — >■ {f,n^). During radiation 

domination, - 
to 
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gcsT"^, Eq- ( 11 1 simplifies further 



dY 

dx 



-\x-'-iY'-Y,l), 



(12) 



where we have assumed Eq. ([t]) to hold and introduced 



A = 



9*s 



1 - 



3 9*s 




mpim^ao ■ 



(13) 



At very early times, Eq. (Ill forces Y to follow the 



equilibrium value Y = Y^q, i.e. WIMPs are very effi- 
ciently kept in chemical equilibrium with the heat bath 
and detailed balance is maintained between DM annihi- 
lation and creation from heat bath particles. At chemical 



decoupling, the factor in front of (1 — Y^^/Y'^) in Eq. pT| ) 



has dropped to a value where DM annihilations are not 
anymore efficient enough to maintain chemical equilib- 
rium and Y starts to decrease at a slower rate than l^q. 
From inspection of this equation, this roughly happens 
when H ^ (crwrci)oq ^^x'^; for WIMPs with the right relic 
density to explain all of the observed DM today, this 
is the case at Xcd — 20 — 28 (see, e.g., [IQ]). Only 
slightly later, the DM density freezes out completely, 
i.e. Y (which roughly corresponds to its comoving num- 
ber density) stays constant - for standard WIMPs essen- 
tially until today. In the case of a Sommerfeld-enhanced 
annihilation rate, the expected range of Xcd essentially 
does not change. Compared to the standard case, how- 
ever, the final stage of the freeze-out process may be de- 
layed. For the models that we will consider here, this 
corresponds to an increase in X99 by up to a factor of 20, 
or even more in the case of resonances, where we define 
X99 as the value of x when Y differs from its asymptotic 
value by less than 1%. 



B. Kinetic decoupling 

After chemical decoupling, WIMPs are still kept in lo- 
cal thermal equilibrium with the heat bath by the much 
more frequent elastic scattering processes with SM par- 
ticles. ^ These contributions to the collision term read: 



1 



d^k 



d^k 



d^p 



2gx^J (27r)32a;7 i2TTf2Cj J {2ttY2E 
x{2nf5('\p + k-p-k)\M\l^^^^ 
lT9^{Lo)) 5^(^)/(p) 



(iTff^('i)) g^Hfiv) 



(14) 



where 4-momenta with (without) a tilde describe ingoing 
(outgoing) particles and \M\^ again is the matrix element 
squared and summed over all spin states. Note that Coi 
only contains particle number-conserving processes, so it 
does not contribute to Eq. Q - as can straight-forwardly 
be checked explicitly [S]. 



^ When including coannihilations this equation takes the same 
form - with Eq. being replaced by an effective thermally av- 
eraged cross-section [cFcBv) and n^'^^ denoting the total (equi- 
librium) number density of all co-annihilating particles [6]. 



^ See Ref. 1361 for an example of how to treat even inelastic scat- 
tering processes in case there exists another heavy particle highly 
degenerate in mass with the DM particle x- 
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Rather than the first moment of the Bohzmann equa- 
tion, as in Eq. Q for the determination of the chemi- 
cal freeze-out temperature, one may consider its second 
moment to get an accurate description of when the DM 
particles leave thermal equilibrium with the heat bath 
[21 [To]. To this end, it is very convenient to introduce 
the parameter 



9x 



(15) 



which would correspond to the temperature of a non- 



relativistiv WIMP if f were a thermal distribution. The 
difference between and T thus indicates how well the 
WIMPs are kept in thermal equilibrium with the heat 
bath. In analogy to Eq. (10 1, we further introduce the 
dimensionless quantity 



y 



- „2/3 



(16) 



Multiplying Eq. ([!]) by g^p'^/E, integrating it over p 
and keeping only the leading order terms in p^/m^ then 
leads, after a somewhat lengthy calculation, to 



with 



(cwrci): 



and dO] 



y_ 
y 



iTm^nl J (27r)3 J (27r)3 



V 3g*s/ Hx \ y J ' 



p' (vroii^^^xx) f{E)f{E) 



SKlix) 



(CTiJroi) VS- 1(2S - 1) [2(5 - l)Ki{2x^/S) + S"5(45 - l)K2{2xVs) 



ds 



^j=- (fJWroi) (^1 + ^OT^ ) e-^'^d?; , 



c(T) 



12(27r 



X X ■'^—rr 



(17) 

(18) 
(19) 
(20) 

(21) 



In arriving at Eq. (17) 



we have assumed that Foq ^ Y 
and Eq. ( 19 ) is valid if f{E) oc e^^l^ (s is the dimension- 
less version of the Mandelstam variable s = 4m^s). The 



approximation given in Eq. ( 20 1 takes a form very simi- 
lar to Eq. ^ and also exhibits an accuracy that is very 
similar, i.e. much better than at the percent-level for the 
values of x that we are interested in here. For a scaling 



like in Eq. ([t]), crvici oc u^", in particular, we find 



( 0-^^1-01)2 

(O-Wrol) 



= 1 



n 
3 



(22) 



Note that, for Y' — 0, Eq. (17) is the exact analogue of 
Eq. ( 11 ): as long as the scattering processes are frequent 
enough, y follows the heat bath value T/eq = 'm^Ts~^^^ , 
i.e. we have = T as expected. At very late times, on 



^ This expression improves the corresponding Eq. (10) in Ref. |10| 
by using a more suitable definition of y; more importantly, we 
allow here explicitly for the case that DM annihilation has not 
ended yet (i.e. dY/dx ^ 0). 



the other hand, the factor in front of (l — j/cq/y) becomes 
vanishingly small and y stays constant, i.e. oc a~^, 
which simply reflects the redshift of the WIMP momenta 
due to the expansion of the universe. The fact that the 
transition between these two regimes happens on a rather 
short timescale [THI allows to conveniently define the 
temperature of kinetic decoupling as 



a^kd 



^2/3 
212 



(23) 



T=Tk, 



As expected, kinetic decoupling happens considerably 
later than chemical decoupling; in the case of neutralino 
DM, e.g., one finds a;kd/a;cd ~ 10 - 4000 (or ^ 
5MeV-5GeV) 10 . 

As we will see, things may change considerably in situ- 
ations where we cannot actually neglect Y' . We therefore 
advocate, as indicated, to use the above definition of a;kd 
only after setting Y' /Y = by hand in Eq. (17 1. This 



definition then accurately reflects the intuitive meaning 
of kinetic decoupling even in the case where we cannot 
neglect Y' , i.e. the point where scattering processes with 
heat bath particles are no longer effective. 
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III. EVOLUTION OF DARK MATTER 
DENSITY AFTER KINETIC DECOUPLING 

In the conventional WIMP scenario, the colhsion term 
m Eq. (g can be completely neglected by the time 
of kinetic decoupling, i.e. the further evolution of / 
is only governed by the expansion of the universe - 
at least until the tiny primordial density fluctuations 
have grown large enough to trigger structure formation 
and self-annihilation may start again. For the case of 
Sommerfeld-enhanced annihilation rates, as we will dis- 
cuss now in some detail, this part of the evolution history 
is qualitatively different and much more complex. 



A new era of annihilation 



to 



T^. For a 



the following Section IIIB), we can use Eqs. depOl ) 
calculate {avrei){2) simply by replacing T ~ ~ 
Sommerfeld enhanced s-wave annihilation, e.g., we then 
have 



ttT. 



f , 



(27) 



where the last step is valid if velocities of the order of 
V ^ V = y/3T^Jrn^ fall into the Coulomb regime where 

S{v) cx v^^; this is exactly the T^^^^ oc x^l'^ scaling men- 
tioned above. For a full understanding of the evolution 
of the WIMP number density and "temperature" in this 
regime, however, we need to solve the following coupled 
system of differential equations for y and Y that follows 
from Eqs. (fTTl[T7|: 



Let us focus on the standard situation where a:kd ^ 
Xcd; around and after kinetic decoupling, we thus have 
Y >• l^q. Therefore, the formal solution to Eq. (Ill is 
given by: 



Y{x)-^ ^Y{xi)-^ + 



3 5*s 



s(crt^rcl) 
Hx 



dx, (24) 



for any x^ ^ x^d- In order to gain some qualitative un- 
derstanding of this expression, let us again assume that 
CTfrci oc w^". As discussed in the previous Section, we 
roughly have v ~ vl'^x °^ before kinetic decou- 

pling and V oc x^^ afterwards; as a consequence, we ex- 
pect (crurci) oc a;~", where 



for X < 



a;kd 



2n for X > a;kd 



(25) 



Approximating A given in Eq. (13) to be constant, we 



can now integrate Eq. (24) and find 



Y(x)-^-Y{xC)-^ 



T^(^-^) forn^-1 
ln(x/a;i) forn = — 1 

(26) 

Clearly, an appreciable change in Y for x > xi \s only 
possible for n < — 1; in fact, taken at face value, annihi- 
lations would never cease in that case. For the standard 
WIMP scenario, this is impossible to achieve since s-wave 
annihilation implies n = h — Q and higher partial waves 
are even more strongly suppressed (e.g. n = \ for the p- 
wave). For a Sommerfeld- like 1/w enhancement of s-wave 
annihilations, however, the situation looks very different 
and WIMPs may re-enter an era of annihilation [17 : in 
this case, we do have h ~ —\ after kinetic decoupling. 
On resonances, we could actually have (crWioi) oc 
i.e. h = —2 (see Appendix [A|) ; note that this would im- 
ply a non-negligible annihilation rate even before kinetic 
decoupling (with n = h — —1). 

Let us now have a more detailed and quantitative look 
at this effect. Assuming that the DM velocity distribu- 
tion stays Maxwellian even after kinetic decoupling (see 



Y' 


1 


3 g,s 


Y 




Hx 


I = 


1 


X g,s 

3 g,s 


y 




Hx 



2m^c{T) 1 - 



y 



(28) 
(29) 



^sY 



This set of equations provides one of our central results; 
it clearly demonstrates that kinetic and chemical decou- 
pling cannot, in general, be treated separately. 

Some insight in the asymptotic behavior of these cou- 
pled equations is achieved by considering the limit where 
X ^ Xkd, i-e. where the scattering term proportional to 
c(T) can be neglected. Assuming again crurci oc w^", and 
using Eq. (22), we then find 



y_ 
y 



nY' 



hY' 
6Y 



(30) 



For n < 0, a decreasing Y will thus have the effect of 
increasing y even after kinetic decoupling; this simply 
refiects the fact that the DM phase-space density is de- 
pleted of low velocity particles, thereby increasing the 
average velocity. 

In the remainder of this Section, we will continue our 
discussion of the further evolution of Y and ?; on a rather 
general level; in Section |IV[ we will then consider a con- 
crete class of WIMP DM models and show that the ef- 
fects discussed here can, indeed, be quantitatively quite 
important in determining the relic density or the small- 
scale cut-off in the mass-distribution of DM subhalos. 



B. Dark matter self-scattering 

In the presence of a Sommerfeld-enhanced annihilation 
rate, al so th e DM self-scattering rate is enhanced (see Ap- 
pendix A 2 1; as a result, WIMPs can have a Maxwellian 



velocity distribution even after kinetic decoupling has 




T [MeV] 
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[TeV] 


[GeV] 


a 


91 


[MeV] 


[MeV] 


off-ros. (min) 


5 


5 


0.0955 


10-^ 


400 


11.5 


off-res. (max) 


0.5 


0.1 


0.0154 


10-1 


0.0665 


4 X lO"'^ 


on-res. (min) 


5 


5 


0.0395 


10-1 


8.87 


Tstruc 


on-res. (max) 


0.1 


0.1 


0.00168 


10"^ 


0.145 


< T' 

^ J-struc 



TABLE I: Model parameters for the minimal and maximal 
cases (with respect to the resulting Tnt) that are shown in 
Fig.[l] For resonances, self-scattering ceases to be effective at 
a temperature Tstruc around the onset of structure formation 
(the effect of which is not included in Fig. [l|. 



FIG. 1: The Hubble rate (full, black) in comparison with the 
effective DM self-scattering rate, introduced in Eq. ( |31[ ), for a 
few extreme parameter sets (minimal in red, maximal in blue) 
of the model considered in Section IV (see also Table |l|. The 
dashed and dotted lines show the on- and off-resonant case, 
respectively. The-dotted line indicates, for comparison, the 
temperature of matter-radiation equality. The temperature 
Tnt below which the DM velocity distribution starts to deviate 
from a Maxwellian form is roughly given by rs(rnt) ~ H{Tnt)- 



taken place [H [371 [3H]- In this already indi- 

cated, we can easily evaluate the thermal averages that 
appear in Eqs. (28 29 1 by using the expressions in Eq. ^ 



and Eq. (20) with T ^ T^. We thus would like to esti- 



mate when the DM self-scattering ceases to be effective 
and the DM velocity distribution starts to deviate from 
a Maxwellian form. 

The self-scattering rate by which the velocities change 
by 0(1) is given by [HEH] 



9l 



rcl 



, , 27r3 2^3 '^TVr.l ^2 



f{E)fiE), 



5/2 



where in the last step / oc e ^/'^x was used. Here, vq = 
y/2T^Jrn^ is the most probable velocity and the transfer 



cross-section ctt is introduced in Eq. (A12| 



Following Ref. [T^, we define the temperature Tnt at 
which the DM velocity distribution becomes non-thermal 
as the temperature where the scattering rate Eg becomes 
comparable to the Hubble expansion rate 



F,(T„t) = Tr(Tnt). 



(32) 



We note that a more precise determination of Tnt would 
in principle be possible by solving the second moment of 
the Boltzmann equation with a collision term that de- 
scribes WIMP self-interactions instead of scattering on 
SM particles. Unfortunately, however, one can no longer 



use the fact that the momentum transfer is small in this 
case and a substantial extension to the formalism pre- 
sented in Refs. [9l[T^ would be required, which is beyond 
the scope of this article. 

For illustration, we show in Fig. [l] the evolution of 
rs(T) and H{T) in the resonant and non-resonant case, 
for two situations that represent rather extreme cases 
for the parameter space of the DM model that we con- 
sider in Section |IV| (see Table |l] for the relevant model 
parameters and the resulting decoupling temperatures). 
As can be seen, the self-scattering rate stays larger than 
the Hubble rate for a very long time if the Sommer- 
feld enhancement is resonant and the scattering rate is 
always able to keep the DM velocity distribution ther- 
mal beyond matter-radiation equality at Tgq ~ 0.75 eV 
[T]. Once matter domination sets in, the Hubble rate 

H (X T'^''^gl's(T) will catch up with F^, but for many 
models it stays below the self-scattering rate all the way 
down to the onset of structure formation (when the for- 
mation of gravitational potentials leads to an increase 
in the DM velocities). If we are not near a resonance, 
the intersection of Eg with H will take place for much 
larger temperatures: in our model, we roughly expect 
0(10) MeV < T„t < O(IOO) eV. Even when following 
a conservative approach for the numerical calculation of 
we find that a thermal velocity 



dv , (31) (see Appendix A 2 



distribution is ensured also for all off-resonance models 
considered here (i.e. before the annihilations finally come 
to an end, see below). 

In summary, we can safely assume a Maxwellian 
WIMP velocity distribution in all cases relevant to our 
discussion - though we stress that this is a nontrivial 
statement that needs to be checked explicitly when ap- 
plying our treatment to other DM models. The most 
noteworthy exception happens if the model parameters 
are tuned in such a way as to very efficiently suppress the 
self-scattering rate due to the Ramsauer-Townsend effect 
[55] . In this case, the analysis of Eqs. ( 28|29 ) is compli- 
cated by the considerably more involved determination 
of the velocity averages that appear in these equations. 
We do not consider this interesting possibility any further 
here but leave it as a challenge for future studies. 



7 



C. Final relic density 

As explained earlier, a new era of annihilations after 
chemical (or even kinetic) decoupling is realized when 
Sommerfeld enhancements are taken into account. We 
now discuss in some detail three effects that can cause 
these annihilations to cease eventually. Obviously, the 
temperature at which this happens is important for the 
correct determination of the relic density. 

The first and most important effect is that the Som- 
merfeld enhancement does not continue to scale like v^^ 
or even as u — > (see Appendix [A|, but saturates 
below some cutoff velocity; at this point, DM annihila- 
tions will no longer be able to keep up with the expan- 
sion of the universe. One can estimate the temperature 
Tsat at which this happens by equating the mean WIMP 
velocity, v = y^ST^/m^ for a Maxwellian distribution, 
to Woff ~ 0.5 m^/m^ (won ^ w^/m^) for a Som- 
merfeld enhancement off-resonance (on-resonance), see 
Eqs. (A8 A9|. Assuming that kinetic decouphng takes 
place more or less instantaneously at rkd7 and that we 
have Tsat < Tkd (which is the situation we are interested 
in here), the saturation temperature can thus be approx- 
imated by: 



(33) 
(34) 



Another effect that is relevant for the discussion of the 
final relic density is the onset of matter domination. We 
can see this directly by using 




47r3 
45to|[ 



<?cffT4 + Tj^Pm,oa-^(r) (35) 



3m|,j 



(36) 



in Eq. (24), where p,„^o = il^Pc 



1.10 X lO^-^^GeV^ is 

the matter density today [1 . Evidently, H^(x) cx x~'^glg^ 

and H^{x) oc x~'^l'^g'Jg , i.e. the effect of matter domi- 
nation on the evolution of the DM density is essentially 
the same as replacing h in {av,-ci) oc x"" by n' = 1/2. 
As described in section |I11 A[ an appreciable change in 
Y is then only possible for n' < —1, which is now no 
longer achieved with a 1/w enhancement after kinetic de- 
coupling (n' = —1/2). Only in the special case that we 
are on a resonance, we would have n' — —3/2 after ki- 
netic decoupling has taken place, and annihilations could 
continue even beyond matter-radiation equality. 

The last effect, which eventually brings all cosmological 
WIMP annihilations to an end, is the onset of structure 
formation around Zgtruc ~ 0(100): once significant grav- 
itational potentials are formed, they start to attract the 
DM particles and cause their average velocity to increase 
again, eliminating the Sommerfeld effect. 



For the models that we are interested in here, to be 
introduced in Section [IV[ we always find 



^ sat, off 



> T„t » T, 



oq ■ 



(37) 



This means that the final end of DM annihilation is trig- 
gered by the saturation of the Sommerfeld effect at small 
velocities as long as the particular combination of DM 
and exchange particle masses do not put us on a reso- 
nance. On a resonance, on the other hand, the saturation 
of the Sommerfeld effect happens much later and we find 
cases in which the WlMPs continue to annihilate until 
well after matter-radiation equality (though we always 

have Tsat, on > Tstruc > Tnt). 

Taking into account the appropriate saturation ef- 
fect (s) described above, let us now denote with Yq the 
result of integrating Eq. (28) to xq — m^/Tp, where 
To = 2.348 X 10~^ eV is the photon temperature of the 
universe today. The DM relic density, in units of the 
critical density, is then obtained as [ID] 

= m^soYo/pc = 2.742 x lO^'h-^ Yo . (38) 

For viable DM models, this should be compared to 



n 



DM 



= 0.229 ±0.015, 



(39) 



which is the most recent estimate from observations [T]. 
If the DM particle is its own antiparticle, X = X: 



should thus demand il-^ 
^^dm/2. 



n 



DM 



otherwise we have fly 



D. The smallest protohalos 

A new era of annihilations after kinetic decoupling not 
only affects the DM relic density, but also the DM veloc- 
ity dispersion at the onset of structure formation. This 
quantity, in turn, can be related to a small-scale cut- 
off Mcut in the power spectrum of matter density fiuc- 
tuations which corresponds to the mass of the smallest 
gravitationally bound objects [Hj . 

Let us therefore use the asymptotic value of ?/, when 
both DM scattering and annihilation events have finally 
come to an end to define an effective asymptotic decou- 
pling temperature 



■^doc 



y\x- 



dec 



5^/3 
2^2 



(40) 



T=T? 



Note that this temperature can be quite a bit higher than 
the temperature at which the annihilations actually stop. 
It thus does not have any obvious intuitive interpretation 
(unless there is no second era of annihilation, in which 
case it is simply given by = Tkd), but is only intro- 
duced here for convenience because it allows to express 
the asymptotic DM temperature, aka velocity dispersion, 

as = ryr^^^. 
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As a consequence, we may simply replace Tkd — >■ 
in any expression that relates the kinetic decoupling tem- 
perature to the cutoff in the power spectrum of matter 
density fluctuations. In particular, we have |10| 



Mcut = max [Mis, Mao] , 



(41) 



where free-streaming 114) induces an exponential cut- 
off characterized by 



Mfe sa 2.9x10" 



1 + ln (seff'l^i^c/SO Mev)/19.1^ 



V 100 GoV 



.1/2 1/4/ Tjg^ \ 
J 'JcS y 50 McV J 



1/2 



(42) 



and acoustic oscillations [151 [IS a similar cutoff with 



1/4^ 



3.4 X 10 



-6 ( -^doci/cff 



50MeV 



Mc- 



(43) 



We note that the above prescription for calculating 
Mcut can, strictly speaking, only be applied if / i) reaches 
its asymptotic behavior (given by y ~ J/L-i-oo) already 
well before the onset of matter domination and ii) has a 
velocity distribution at that time which resembles that of 
the standard WIMP case - simply because these are the 
assumptions that entered into the analysis of Refs. [TTl 
[TlHIl]- As we will see, it is possible to find viable models 
where these assumptions are not satisfied; in these cases, 
our prescription thus only gives an approximative value 
for Mcut- While we do not pursue this issue any further 
here, it would certainly make for an interesting study 
to extend the standard analysis and investigate how the 
spectrum of matter density fluctuations would be affected 
in such cases. 



IV. LEPTOPHILIC DM MODELS 

In this Section, we demonstrate how to apply our gen- 
eral discussion to a specific class of WIMP models that 
exhibit large Sommerfeld effects. As already mentioned 
in the introduction, the cosmic ray lepton excess 31j has 
triggered a lot of phenomenological activity, trying to es- 
tablish a possible DM connection. At first sight, this 
is not a trivial task at all since standard WIMP can- 
didates, like the neutralino, fail to meet the necessary 
criteria from a model-independent analysis of the data 
[411 l42l in order to fit the observed excess, the annihila- 
tion of DM particles would have to i) produce a very hard 
positron spectrum at least up to a few hundred GeV, but 

ii) at the same time produce almost no antiprotons, and 

iii) happen with a rate about 2 to 3 orders of magnitude 
above the canonical value of 3 • 10~^^cm^s~^ which is 
needed for thermal production in the simplest scenarios 
(see, however, Ref. [43] for a way of how neutralino DM 
could at least satisfy the two first requirements). 

An elegant, or at least very economic, way to meet 
all these requirements is to postulate the existence of 



a new light exchange boson (j>, with 100 MeV < < 
1 GeV, that directly couples to a DM particle with ~ 
0(1 TeV) but only very weakly to standard model par- 
ticles (in order to avoid the stringent bounds on light 
new particles). The DM annihilation rate today would 
then be strongly enhanced by the Sommerfeld effect, 
without changing too much the annihilation rate dur- 
ing freeze-out; furthermore, the decay of the resulting 
particles into hadronic modes is kinematically forbidden, 
thus leaving only the desired light leptonic modes (see 
Ref. [53] for a general account of this idea). 

Many realizations and variations of this rather general 
setup have been worked out in, e.g., Refs. [2SH2^ . Noting 
that an astrophysical explanation of the above mentioned 
cosmic ray excess may well be more likely |44] . we will 
here only consider a very simple phenomenological toy 
model for leptophilic WIMP DM - loosely motivated by 
the cosmic ray lepton data but mostly chosen for illus- 
tration of the possible effects of the Sommerfeld enhance- 
ment on the thermal decoupling process. 



A. A simple toy model 

We consider a fermionic DM particle x that couples 
only to a light scalar (j)s and a pseudo-scalar (j)p via 



(44) 



We assume that the (pseudo-)scalar particles also inter- 
act with standard model leptons, through 



(45) 



g^'f = and 



albeit with much smaller coupling strengths (i.e. and 
gf can be thought of as effective couplings arising, e.g., 
from higher-dimension operators). 

For simplicity, we will assume gl^^ 
gl'P = 0; possible couplings to quarks typically do not 
change the phenomenology of our model (for the mass 
ranges considered here) and are anyway strongly con- 
strained - see Appendix |B] for experimental bounds on 
this and similar models. The scalar particle mediates DM 
scattering at low momentum transfer and is also respon- 
sible for the Sommerfeld enhancement. The need for an 
additional pseudoscalar particle arises because of parity 
conservation: While the annihilation XX ~^ 4's4'p has an 
s-wave contribution, xX ~^ 4>s4's{4>p4>p) is a p-wave pro- 
cess which vanishes in the w — ;> limit, so it is only the 
former channel that can happen with a sizable rate today 
(which, in turn, is required if we want to make any con- 
tact to observable cosmic ray lepton fluxes at all). Even 
the relic density is to a large extent determined through 
this channel and thus mainly depends on the parameter 



An 



(46) 



Since we want to keep the discussion simple, we will 
mostly assume 5* = 5^ in the following, noting that 
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FIG. 2: Shown are two approximations (dotted, dashed) 
to the numerically calculated, full Sommerfeld enhancement 
(solid), for = 1 TeV, = 5 GeV and a = 0.03. See text 
for further details. 



the relic density is essentially only set by one effective 
coupling constant anyway. 

We note that our setup roughly corresponds to a sim- 
ple version of the model proposed by Nomura and Thaler 
[28j . where the new particles postulated above (x, 0s, 4>p) 
are embedded in a full supersymmetric scenario and 
4>p, in particular, takes the role of an almost 'standard' 
Peccei-Quinn axion. In Appendix [C] we collect the rel- 
evant annihilation cross-sections and scattering matrix 
elements for this model. 



B. Decoupling on and off resonance 

For this model, let us first study in some detail the 
evolution of the WIMP number density and temperature 
for two cases of particular interest, i.e. parameter sets for 
which we are outside and exactly on a resonance, respec- 
tively. We will be especially interested in quantifying the 
difference between our full treatment, described in Sec- 
tion[II|and |III[ and various approximations one may deem 
reasonable. In particular, we will refer to the full solu- 
tion as the coupled set of Eqs. (11 17), with the proper 
replacement of T — Ty in the expressions for (trwr, 



for Y ^ yoq, this simply corresponds to Eqs. (28 29 1 



1 /(2 



in Eq. (11 



Eq. (171 



We denote with sudden decoupling the case in which we 
assume 7^(T > Tkd) = T and T^{T < Tkd) = rVTkd 
) and neglect the term proportional to Y' in 
?he uncoupled solution, finally, corresponds to 
the case where kinetic decoupling is assumed to have no 
influence on the evolution of the WIMP number density, 
and vice versa. 

In addition to the full numerical solution for the Som- 
merfeld factor, we also considered two different approx- 
imations as illustrated in Fig. [2] For approximation 
1 (shown as a dotted line) we assumed S = 1 for 
V > na, S = na/v in the intermediate regime and 



S{v < Vmax) = S'max (defined by the full solution). As a 
better approximation, especially at larger velocities, we 
used the Coulomb approximation S = ^ /{I — e~^) 
down to Wmax, and S — S'max for smaller velocities. This 
approximation 2 is shown as a dashed line. For the case 
of resonances, we adjusted these approximations corre- 
spondingly: in approximation 1 we used S cx v~'^ instead 
of 5* cx v^^, and in approximation 2 we used the Coulomb 
expression down to velocities where the enhancement fol- 
lows the behavior. 

In Fig. |3j the solutions for y and Y are shown for one 
particular parameter set for which the Sommerfeld en- 
hancement is not in the neighborhood of a resonance. 
As expected, the relic density is lower (by 10%) for the 
coupled solution than for the uncoupled solution; com- 
pared to the case without Sommerfeld enhancement, it 
is smaller by 80%. Also visible is that approximation 
2 reproduces the full numerical result much better than 
approximation 1, for which the relic density is overesti- 
mated. Assuming that kinetic decoupling happens sud- 
denly, like depicted in the right panel of Fig. |3] in red, 
gives a reasonable approximation to the full evolution 
of y - though we find that for very small lepton cou- 
plings, and negligible (/)-scattering, the kinetic decoupling 
process may be more delayed than for standard WIMPs 
[5J[TD]: since more scattering events are necessary to keep 
the DM in thermal equilibrium in this case, it also takes 
longer before scattering becomes sufficiently inefficient so 
that y remains completely constant. Therefore, sudden 
decoupling works very well as an approximation to Y 
(the lines are completely obscured by the full solution 
in Fig. [sj right). For this set of parameters, kinetic de- 
coupling happens at ^ 2.92 x 10"^, after which the 
coupled solutions start to deviate from the uncoupled 
solutions. The annihilations cease around a; 3 x 10^ 
because of the saturation of the Sommerfeld factor, and 
both y and Y stay constant afterwards. 

The more interesting case is to consider a parameter 
set for which we are on a resonance in the annihilation 
cross-section. The results, shown in Fig. |4j show a strik- 
ing difference with respect to the uncoupled Boltzmann 
equations for x > 10^, after which the WIMP annihila- 
tions decrease the relic density by more than two orders of 
magnitude until the Sommerfeld enhancement saturates 
at X ~ 2 X 10^'^. Following the argument presented in Sec- 
tion |III K\ one might have expected an even more efficient 
decrease in the relic density that starts directly after ki- 
netic decoupling (which for this choice of parameters hap- 
pens at — 2.35 X lO**). This is not observed, however, 
because at early times the velocities of the WIMPs are 
still large enough to be in the Coulomb regime and, even 
though we are on a resonance, the Sommerfeld enhance- 
ment follows an S' cx 1/v rather than S ex. 1 jv^ behavior 
(see Appendi x [A| . Indeed, for 3 x 10^ < x < 3 x lO*" 



find that Eq. (30) is satisfied with — 2 < n < — 1, whereas 



the Sommerfeld enhancement shows the resonant behav- 
ior (n — —2) for larger values of x. 

The value for the coupling ol is here chosen such that 
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FIG. 3: For a parameter set where the Sommerfeld enhancement is not near a resonance {rrix = 1 TeV, = 5 GeV, a = 0.03 
and gi = 10~^), we show the evolution of the quantities Y as defined in Eq. ( |10[ |, and y as defined in Eq. ( |16| |. We show the 
solution to the full set of coupled Boltzmann equations (black) as well as for two approximations described in more detail in 
the text, assuming sudden kinetic decoupling (red; not visible in left plot) and no coupling between the Boltzmann equations 
for y and Y (blue) , respectively. Different approximations to the Sommerfeld enhancement are shown according to Fig. [2] by 
dotted, dashed, or solid lines. Additionally, the solution without Sommerfeld enhancement is shown in the left panel (green, 
dash-dotted). 
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FIG. 4: Same as Fig. [3] but now for a parameter set where the Sommerfeld enhancement is resonant (m^ = 1 TeV, = 1 
GeV, a = 0.00168 and gi = 4.6 x lO"''). It can be seen that annihilations continue even until after matter-radiation equality 
(denoted as a dash-dotted line). 



the full solution gives the right relic density today (within 
3(t): the uncoupled solution (both with and without 
including the full Sommerfeld factor) actually overesti- 
mates Yb by a factor of ^ 400. Sudden decoupling be- 
comes a very bad approximation to y in the case of res- 
onances, at least for x > 10*" where it simply follows the 
uncoupled solution. However, it catches the overall be- 
havior of a large decrease in Y rather well - though it 
underestimates the final relic density by at least a factor 
of 2. The different approximations to the Sommerfeld en- 
hancement, on the other hand, give comparable results. 



C. Range of decoupling temperatures and the mass 
of the smallest protohalos 

Having discussed in some detail the situation for two 
particular parameter sets, let us now explore the ther- 
mal history of DM, and in particular the consequences 
of our improved treatment, for the full possible range of 
our model parameters. For this purpose, we adjusted 
the coupling a in all calculations in such a way that the 
relic density obtained by solving the full, coupled Boltz- 
mann equations is within Scr of the observed value today, 
0.184 < r^DM < 0.274, leaving thus the lepton coupling 
gi as the only other free parameter besides and m^. 

Let us start by showing in Fig.[5]the kinetic decoupling 
temperature T^d as a function of the coupling constant ge- 
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FIG. 5; For models where the Sommerfeld enhancement is not 
on a resonance, the kinetic decoupling temperature is shown 
as a function of the mediator particle coupling to leptons, for 

= 100 MeV (black, fuU), 500 MeV (blue, dashed), 1 GeV 
(red, dotted), and 5 GeV (green, dash-dotted). From bottom 
to top, the lines correspond in each case to a DM mass of 

= 100, 500, 1000, 5000 MeV. 



FIG. 6: For models where the Sommerfeld enhancement is 
resonant, the kinetic decoupling temperature Tkd (filled) in 
comparison to the effective asymptotic decoupling tempera- 
ture T^;. (empty) is shown. Again — 100 MeV (black, 
squares), 500 MeV (blue, circles) and 5 GeV (green, trian- 
gles). 



As can be seen, Tkd decreases for larger gi - reflecting the 
fact that a strong lepton coupling will keep the WIMPs 
longer in local thermal equilibrium. For decoupling tem- 
peratures smaller than around 7 MeV, only DM scatter- 
ing with electrons is effective just before the DM particles 
completely leave thermal equilibrium. At higher temper- 
atures, also muons start to contribute very efficiently to 
the scattering process, resulting in a flattening of Tkd(g£) 
when moving to smaller values of gi; this already hap- 
pens for highly non-relativistic muons because of the mf 
dependence of the scattering rate, see Eq. ( |C4[ ). Scatter- 
ing with mediator particles ensures that the kinetic de- 
coupling temperature does not increase arbitrarily high 
even for negligible lepton couplings, which explains the 
plateau that appears at roughly gi < 10~^. The de- 
coupling temperature Tkd increases rather strongly with 
higher mediator masses m^: in the case of lepton scat- 
tering this directly follows from the form of the scatter- 
ing matrix element, Eq. (C4), while in the case of scat- 



tering with the mediator particles it reflects the strong 
Boltzmann suppression of the latter. The dependence of 
the decoupling temperature on the DM mass , on the 
other hand, is very weak. Also the spread in Tkd (for 
given values of m^,m^ and gi) due to the different val- 
ues of the relic density that were obtained by changing 
the coupling a accordingly is essentially negligible. 

In Fig. [5] we have chosen to only show models with pa- 
rameter sets where the Sommerfeld enhancement is not 
resonant both because this is in some sense the generic 
behavior of the model and because the parameter depen- 
dence is very straight-forward to discuss in that case. 
Phenomenologically, on the other hand, the resonant 
case is also very interesting, not the least because much 
stronger effects compared to the standard scenario can 



be expected. We therefore performed a dedicated scan 
to sample that part of the parameter space that results 
in a resonant Sommerfeld enhancement and the correct 
relic density (in practice, we only sampled the first 5 reso- 
nances; higher resonances are more and more densely dis- 
tributed, increasing thus the amount of parameter fine- 
tuning that is necessary to meet the just mentioned crite- 
ria). As shown in Fig.[6j the resulting kinetic decoupling 
temperature is a bit higher than in the off-resonance case. 
This is due to the fact that the relic density constraint 
makes a smaller coupling a necessary to compensate for 
the large Sommerfeld effect, which in turn decreases the 
amount of DM scattering off ips so that WIMPs decouple 
slightly earlier than in the off resonance case. The asymp- 
totic decoupling temperature T^^, on the other hand, is 
up to a factor of ~ 5 smaller than Tkd and always smaller 
than in the off-resonant case (note that off resonances 
and Tkd differ, as expected, by only 3% at most). 
From the asymptotic decoupling temperature, we also 
calculated the corresponding mass of the smallest grav- 
itationally bound objects; the off-resonance results are 
shown in Fig. [7] It can be seen that the possible cut- 
off mass spans a wide range of Mcut/MQ ~ O(10~^° 
- 10'^). Even when taking into account existing con- 
straints on the lepton coupling, see Appendix |Bj much 
larger cutoff masses than in the standard WIMP case 
thus seem possible; imposing gi < 10~^ e.g. results in 
-^cut ^ Mq (whereas A/cut ^ 10"'^ M© for neutralino 
DM fTO^). Compared to Tkd and T^^, M^nt can have a 
stronger dependence on m^, visible by the spread of the 
bands for gg < 1Q~'^. This is due to the fact that in this 
region of parameter space the free-streaming mass Mfs 
(that depends, among other things, on m^), dominates 
over the acoustic oscillation cutoff mass Mao- In the reso- 
nant case, the resulting Afcut is in general higher; naively 
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10"^ 10 - 10" 



FIG. 7: For models where the Sommerfeld enhancement is 
not on a resonance, the cutoff mass is shown as a function 
of the lepton coupling, for = 100 MeV (black, full), 500 
MeV (blue, dashed), 1 GeV (red, dotted), and 5 GeV (green, 
dash-dotted). From top to bottom, the lines correspond in 
each case to a DM mass of = 100, 500, 1000, 5000 MeV. 
On resonance, the value of Mcut can be larger by a factor of 
up to 2 (20) for very large (small) values of g^. 



using Tkd instead of to calculate it would result in 
values up to 2 orders of magnitude smaller. As expected 
from Fig. |6j the effect of resonances is most pronounced 
for small values of Mcut7 corresponding to large values 
of Tkci; as a result, the lowest possible cutoff mass shifts 
from - 3 X lO-^^M© to - 7 x lO^^M©, while the largest 
possible value shifts from 6OOM0 to ~ IIOOMq. 

Let us now discuss the evolution of the WIMP number 
density and the resulting relic density. Annihilations off 
resonance finally come to an end for temperatures rang- 
ing from Tsat ^100 MeV to ~ 1 keV, where Tgat was de- 
fined by F(Tsat)/^o = 0.99. Whereas the maximum satu- 
ration temperature is more or less independent of m^, the 
minimum Tgat is significantly smaller for small mediator 
masses. This comes from the fact that the Sommerfeld 
enhancement increases for smaller values of m^/(m^a), 
so the new era of annihilations after kinetic decoupling 
can last longer for smaller m^; furthermore, kinetic de- 
coupling happens later in this case. We also find that, 
as expected, Tgat decreases with decreasing gi, i.e. with 
increasing Tj^-d. Finally, we observe that for many models 
with large lepton couplings we have Tgat > Tkd, indicat- 
ing that the effect of reentering an era of annihilations is 
negligible, whereas for small gi we find that T^d/rgat can 
reach values up to 0(10"'). 

For non- resonant Sommerfeld enhancement, the re- 
sulting effect of the new era of annihilations on Y is 
significant but never extremely large: the ratio of the 
relic density obtained by the uncoupled Boltzmann equa- 
tions with that obtained by the full solution is at most 
f^x,u/^x,c ^1-1 for the range of parameters that we have 
scanned. This ratio tends to be slightly larger for higher 
WIMP masses, refiecting the fact that a larger Sommer- 





off resonance 


on resonance 




~ 24 - 27 


~ 20 - 25 


Tkd [MeV] 


~ 0.07 - 400 


~ 0.09 - 450 


Td^c [MeV] 


~ 0.07 - 400 


~ 0.06 - 170 


Teat [keV] 


~ 1 - lO'^ 


~ 10"'' - 10"^ 




~ 1 - 1.1 


~ 3.5 - 670 


Mcut[M0] 


~ 3 X 10"^° - 600 


~ 7 X 10"® - 1100 



TABLE II: Overview of the various decoupling temperatures 
in our toy model, for the parameter range that we considered 
here. Also stated is the resulting change in the relic density 
when fully taking into account the coupled Boltzmann equa- 
tions Eqs. (11 [17 1 rather than ignoring the possible impact of 
kinetic decoupling; finally, we summarize the possible range 
of cutoff masses (see also text and Fig. [7|. 



feld effect due to an increased value of a is expected in 
that case (recall that the relic density is set by the an- 
nihilation rate which scales, roughly, as {av) cx a^/m^). 
We also find that, as expected, fi^j-^u/ri^^c grows with in- 
creasing Tkd: the new era of annihilations simply lasts 
longer when kinetic decoupling takes place early. 

In the resonant case, the change in Y by considering 
the coupled Boltzmann equations rather than the un- 
coupled equation is much more significant, yielding a 
relic density that can be up to a factor of ~670 smaller. 
In most cases annihilations even continue after matter- 
radiation-equality, further decreasing the relic density in 
some extreme cases by up to a factor of 4 before reach- 
ing its final value: r2^(Toq) < 4f2^(To). The relic den- 
sity correspondingly saturates at very low temperatures, 
0(10) eV > Tsat > O(10"3) eV. Since the annihilations 
are more efficiently suppressed after the onset of matter 
domination, the relic density could probably not continue 
to decrease a lot more, also because of the already men- 
tioned increase in the WIMP velocity once gravitational 
potentials build up and cosmological structure formation 
starts. The effect we see here due to resonances is there- 
fore (close to) the maximum effect on the relic density 
that can be expected, and increasing the Sommerfeld en- 
hancement (e.g. by an even denser sampling of the model 
parameters near the resonances) are unlikely to have a 
large impact. 

For convenience, we summarize in Tab. [ll] the above 
discussion by showing the range of possible decoupling 
temperatures that we have encountered in our scan of 
model parameters, along with the resulting change in the 
relic density and possible values for the mass of the small- 
est protohalos. 



V. DISCUSSION 

In this Section, we return to some technical issues that 
we have not addressed explicitly in (sub-) Section, indi- 
cate possible extensions to our analysis and discuss pos- 
sible consequences of our findings. 
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One such comment concerns the scattering of DM with 
mediator particles which we have tacitly assumed to be 
in thermal equlibrium in the above analysis: we actually 
checked that this is always satisfied for T > Tkd because 
of the very efficient (inverse) decay processes O £~^£~ 
(for more details, see Appendix [C]) . If one would neglect 
X(f> ^ x4> scatterings, and only take into account DM 
scattering with standard model leptons, the kinetic de- 
coupling temperature would continue to increase with- 
out any bounds for smaller values of gf . As an interest- 
ing consequence, it is thus only due to the presence of a 
thermal population of particles that a situation with 
^kd > Ted cannot occur (note also that the thermal pro- 
duction mechanism for x actually requires the exchange 
particles to be in thermal equilibrium around Ted! this 
implies a lower bound on gf such that one essentially 
cannot get around this argument even for a much wider 
class of DM models than considered here) . 

Let us also mention that in Section [IV[ for simplicity, 
we only took into account DM scattering with scalar par- 
ticles. Including also the scattering with pseudoscalars 
affects only the kinetic decoupling temperature directly: 
other observables, such as the relic density, change indi- 
rectly because of the coupling between Y and y. In this 
case, we find that Tkd for g^ < 10^"* becomes as expected 
a bit lower (at most ~ 20% for gf — g^ and rris — mp), 
with a correspondingly smaller impact on, e.g., the relic 
density. 

Furthermore, we have in our discussion always as- 
sumed that the DM particles couple with the same 
strength to both scalars and pseudo-scalars, i.e. — 
g^. In principle, however, the coupling strength a, 
cf. Eq. (46), that roughly sets the relic density is not 

'^/47r that gov- 



9; 



the same as the coupling strength 
erns the Sommerfeld enhancement. Allowing for 7^ gP. 
would thus imply a larger range of possible a* that is 
consistent with the relic density requirement. While we 
choose not to explore the full phenomenology of this 
option here, we note that the implications for the off- 
resonance case are not expected to be sizable because a 
relatively large value of a severely restricts the range of 
a* (neither nor gP, must be too large in order to re- 
main in the perturbative regime) ; if g^ and g^. are varied 
independently, on the other hand, less fine-tuning is in 
some sense required to arrange for a resonant Sommer- 
feld enhancement. In order to make the first remark a bit 
more quantitative, let us consider the situation of keep- 
ing a fixed, but choosing g^ — 2gP {gP = 2gJ): for the 
same range as considered before for all other parameters, 
the relic density decreases in this case by at most ^-^20% 
(^-^2 %) and the kinetic decoupling temperature increases 
(decreases) by up to ~3% (~10%). 

Allowing for a non-negligible coupling of the exchange 
particles to quarks, g^ ^ 0, would introduce the pos- 
sibility of further DM scattering processes XQ ^ XQ 
at high temperatures. However, since kinetic decou- 
pling usually happens after the SU(3) phase transition 
at Tqcd ~ 170MeV |45], there are no free quarks around 



anymore and this does not affect our determination of T^d 
- unless one adopts rather large values for > 5 GeV 
(if kinetic decouplings happens at T < Tqcd, the cal- 
culation of Tkd is also affected by the dependence of the 
effective number of relativistic degrees of freedom during 
the transition; in our calculations, we used the values 
provided in Ref. |46j). 

We caution again that our results should be interpreted 
with care in some special regions of the parameter space 
where, due to the Ramsauer-Townsend effect, the as- 
sumption of a Maxwellian velocity distribution after ki- 
netic decoupling may not be a good approximation any- 
more (see the discussion at the end of Section III B ) . The 



same goes for the calculation of Mcut in those strongly 
resonant cases where the asymptotic ('free') behavior of 
the WIMP phase space distribution is only reached af- 
ter matter-radiation equality; cf. the disclaimer at the 
end of Section IIII Dl The fact that we find viable models 
where this indeed happens provides a strong motivation 
for future analyses to actually establish an exact relation 
between decoupling temperature and cutoff mass even 
during and slightly before matter domination - corre- 
sponding to Eqs. ( 42|43 ) which, strictly speaking, only 
hold approximately in this case. 

Let us finally briefly discuss possible ways to directly 
probe the cutoff mass - which would provide a fascinat- 
ing new window into the particle nature of DM. While 
gamma rays (mostly via XX ~^ i~^i~"f) from individual 
subhalos close to Mcut are unlikely to be resolved [47] . 
the one-point probability function of the diffuse gamma- 
ray flux [?5| could provide a better future probe; for very 
large cutoff masses, Mcut ^ M0, also other anisotropy 
probes could be sensitive enough [49] . The last comment 
also holds for probing the smallest halos with gravita- 
tional lensing, especially for future astrometric microlens- 
ing missions with unprecedented sensitivities ,50] : in fact, 
it has been argued that even sub-solar objects could cre- 
ate observable strong gravitational lensing events, espe- 
cially when making use of multiple images of time- varying 
sources [5T] (though this will be especially challenging in 
the case we are interested in here because the Einstein 
radius of DM subhalos is much smaller than their virial 
radius) . 

The maybe most interesting indirect probe of the cutoff 
scale, on the other hand, currently comes from the obser- 
vation that galaxy clusters maximize the enhancement of 
a generic DM annihilation signal due to the presence of 
substructure [32l |33] . Assuming that one can extrapolate 
the results of numerical TV-body simulations of gravita- 
tional clustering from the current resolution limit down 
to the much smaller cutoff values in the power spectrum 
that we are interested in here, it was shown in these refer- 
ences that small values of Mcut are strongly constrained 
by gamma-ray observations of galaxy clusters in models 
where the DM annihilation is Sommerfeld-enhanced; in 
the case of leptophilic models, the main source of gamma 
rays would be inverse Compton scattering of the high- 
energy leptons from xX ~^ off cosmic microwave or 



14 



starlight photons. For a very specific model, with param- 
eters chosen such as to result in a positron spectrum that 
could account for the cosmic ray excess (m^ = 1.6 TeV 
and (j) — >■ /i^, e^, tt^ with a ratio of | : | : the authors 
could show that for M^^t = lO^^Mg (IO^'A/q) the maxi- 
mally allowed Sommerfeld enhancement would be 5 ^ 5 
(800).^ Here, S = S{v ~ 124km/s) and it was assumed 
that the annihilation cross section with Sommerfeld en- 
hancement is given by crw = 3 • lO^^^cm'^/s. 

While a detailed analysis that extends the above re- 
sults to more general models with Sommerfeld enhance- 
ments is certainly warranting, it is beyond the scope 
of this work. Here, we simply point out that we have 
presented a way to accurately calculate the cutoff mass 
in any given such model; limits on Mcut can thus di- 
rectly be translated into limits on the model parameters 
(which determine the size of S) . What is most interesting 
in this respect is that these limits will be strongest for 
small values of the lepton coupling gf because this leads 
to smaller values for Mcut - which means that gamma- 
ray constraints from galaxy clusters probe the parameter 
space of leptophilic models from a completely different 
direction than other experiments (see Appendix [b] for a 
comparison). 



VI. SUMMARY AND CONCLUSIONS 

As already extensively discussed in the literature, the 
impact of Sommerfeld-enhanced annihilation rates on the 
DM relic density can be sizable both before [5^H57] and 
after [TTHTO] kinetic decoupling. In this article, we have 
introduced a general framework that, for the first time, 
allows to consistently describe situations where DM an- 
nihilation continues after chemical decoupling and in- 
terferes with kinetic decoupling, improving thus on the 
general praxis of using the kinetic decoupling tempera- 
ture as an essentially free parameter. The coupled set of 
Eqs. 



(28 



29 ) that describes the evolution of the WIMP 



number density and temperature thus provides one of 
the central results of this article. As a consequence of 
our discussion, we have also refined the usual definition 
of kinetic decoupling in order to discriminate it more 
clearly from the point where WIMP (self-)interactions 
finally come to an end, which may happen much later. 

Applying our formalism to a simple leptophilic toy 
model, we find that the impact of a new era of DM anni- 



* This would lead to the conclusion that this or very similar lep- 
tophilic models cannot account for the cosmic ray lepton ex- 
cess without violating the gamma-ray constraints from clusters 
if Mcut < 10*Mq 1331 . While it is certainly not the prime pur- 
pose of our article to make such a connection, let us just point 
out that we could accommodate larger values of Mcut even in 
our simple toy model by choosing tUb <SC 100 MeV; a possibly de- 
sired decay into /x^ and tt^ could in that case exclusively happen 
through the (correspondingly heavier) pseudoscalar tpp- 



hilation (i.e. after kinetic decoupling) on the relic density 
can be significant. Ojf resonance, this effect is at most 
~ 10% (at least in our model) and thus smaller than 
claimed in, e.g., Ref. [T7]. On resonance, on the other 
hand, we have demonstrated that DM annihilation can 
continue until well after matter radiation equality, deplet- 
ing the DM abundance by more than two orders of magni- 
tude after kinetic decoupling. This is a rather new result 
which completely changes the naive picture of associating 
the relic density of thermally produced DM to processes 
restricted to temperatures around Ted ~ w^/25, i.e. the 
very early universe. 

Concerning the cutoff in the power spectrum of matter 
density fluctuations, we find that the resulting smallest 
DM protohalos form with masses in the range of roughly 
Mcut O(10-i"A/q) - 0(10 M©), depending on which 
experimental limits (see Appendix Q and which model 
parameters one chooses to adopt, ^s it turns out, the 
correctly determined cutoff mass can be almost two or- 
ders of magnitude larger than in the case where the im- 
pact of DM annihilation after kinetic decoupling on the 
evolution of the DM phase-space distribution is not taken 
into account. In general, much smaller kinetic decoupling 
temperatures, and thus larger cutoff masses, are possible 
than for typical WIMPs like neutralino or Kaluza-Klein 
DM pUl [55] , which may eventually even help to distin- 
guish between these types of DM models (note that even 
Mcut ^ 10 M0 is possible for exchange particles lighter 
than 100 MeV). Existing limits from gamma-ray obser- 
vations of galaxy clusters _33J may already now be used 
to rule out the smallest values for the cutoff mass; since 
early kinetic decoupling happens for small lepton cou- 
plings, this places limits on the parameter space that is 
complementary to bounds from g — 2 measurements or 
beam dump experiments. 

While we have chosen a specific, and rather simple, 
leptophilic toy model for illustration, let us stress once 
again that our formalism can be used for any model 
with annihilation rates that are enhanced at small ve- 
locities - including more realistic models motivated by 
the cosmic ray anomalies (see e.g. Ref. 59^ for a recent 
discussion), more classical WIMPs like heavy neutralino 
DM |52l [56l [57] or models where the annihilation rate 
is enhanced through the formation of bound states [60] 
or an s-channel resonance |61j . In fact, in order to ob- 
tain reliable estimates for both the relic density and the 
small-scale cutoff of matter density perturbations in these 
cases, we have demonstrated here that it is mandatory to 
use a framework that consistently takes into account the 
intertwined nature of WIMP annihilation into, and scat- 
tering with, heat bath particles in the early universe. 
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Appendix A: Sommerfeld enhancement 

The Sommerfeld effect [20] arises when non-relativistic 
DM particles interact via the exchange of force carriers 
that are much lighter than the DM particles themselves, 
^ m-^. Multiple (^-exchanges then result in non- 
perturbative corrections that can enhance both the DM 
annihilation (see, e.g., pTH5D] ) and self-scattering P^[571 
155] cross-section significantly. An effective resummation 
over all corresponding ladder diagrams is performed by 
solving the Schodinger equation for the two-body system. 



1 

2fi 



- V{r) V, 



where ipk is the full two-body wave-function, k — m^v 
is the momentum of each particle in the center of mass 
frame, /x = is the reduced mass of the system and 

the potential V{r) depends on the relative distance be- 
tween the two particles r. 



1. Annihilation 

Because the potential typically has a much longer 
range than the very short relative distance r at which the 
annihilation takes place, these two effects do not inter- 
fere: the only relevant effect of the potential is to change 
the free wave-function at r = 0. The full annihilation rate 
is thus simply obtained by multiplying the unperturbed 
rate with the enhancement factor 



siv) = \M0)f 



(A2) 



where 4'k is the solution to Eq. (]A1|) and normalized as 

OD . (A3) 



■0 — >■ exp 



ikz 



exp 



\kr 



for 



For a massless force carrier the potential is Coulomb- 
like, V{r) — —a/r, and the Schrodinger equation can be 
solved analytically to give 



S{v) 



7r/e„ 



1 



-■K /e^j 



(A4) 



where = v/a. For small velocities, e^, ^ 1, the Som- 
merfeld enhancement thus is given by 5 ~ tt/eu, while 
there will be no enhancement for large velocities, i.e. 
S" sa 1 for e„ » 1. 

When the mediator has a non-zero mass m^, the 
potential becomes Yukawa-like, V{r) = —a/r e"'"*'", 
and Eq. (Al) can only be solved numerically (see, e.g.. 



Ref. [55] for a pedagogic treatment). Since V{r) is ro- 
tationally symmetric, we can expand '0fc into products 
of Legendre polynomials and radial functions R^i. The 
only Rki that is non-zero at r = has angular momen- 
tum I — 0, and therefore we substitute Rko = £,/r. ^ then 
obeys the radial Schrodinger equation 



with the following boundary conditions: 

m - 0, 



dr 



= im-^v ^ for r 



(A5) 

(A6) 
(A7) 



Expanding the Yukawa potential in powers of m^r, 
V ~ —a/r + am^ + ©(r^), one recovers the S chro dinger 
equation with a Coulomb potential for Eq. ( A5 1 when 



(Al) am(f) <^ m^v . Therefore one can safely use Eq. (A4) as 



a good approximation to S for a massive mediator in the 
range > = m^/{am-^). 

For smaller velocities, the Coulomb approximation is 
no longer valid, i.e. the exchange particle can no longer 
be treated as effectively massless. In this case, the fi- 
nite range of the Yukawa potential induces resonances 
in S that correspond to quasi-bound states of the pair 
of the incoming DM particles. We can get some more 
insight in the behavior of these resonances by consider- 
ing the Hulthcn potential instead, which can be tuned to 
approximate the Yukawa potential very well and has the 
advantage to be analytically solvable [53]. In this case, 
the resonances appear at specific values of ~ 6(7m)~^, 
where rt is a positive integer, and the Sommerfeld en- 
hancement at a resonance is given by S* ~ (7r^/6)e0/e^ 
[T^ IM] . In reality, the Sommerfeld factor is limited by 
the finite lifetime of the bound state and therefore should 
not diverge as w — 0; rather, it is expected to saturate 

for [niini 



(A8) 



In our treatment 



we therefore always use S{ev < 

Finally, we note that off the resonances, the full nu- 
merical solution to S saturates much earlier and stays 
constant for |30j 



(A9) 



This simply reflects the fact that the total effective energy 
seen by a Coulomb-like potential, Ecs ^ am^ + m^v^ 
from the expansion of the Yukawa potential, essentially 
no longer depends on the velocity in this regime. 



2. Self-scattering 

In the case of DM self-scattering, the repeated ex- 
change of a gauge boson does not only affect the initial. 
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but also the final state wavefunction; a separation be- 
tween short and long distance scales is thus no longer 
possible and, in principle, all angular momentum con- 
tributions need to be taken into account. The radial 
Schrodinger equation describing self-scattering is thus the 
same as in the case of annihilation, Eq. (A5), but with 



an additional centrifugal term 1(1 + l)/(m^r^) in the po- 
tential. Its numerical solution can be found by imposing 
appropriate boundary conditions, i.e. demanding that 
is regular at the origin and 



4/ sm(m^vr - y + 0;) 



for 



(AlO) 



From the phaseshift Si , one can then calculate the partial 
scattering cross-section as 



47r 



{21 + 1) sin^ Si . 



(All) 



For our purposes, we are mostly interested in the trans- 
fer cross-section 



aT = dfl {1 — cos ( 



(A12) 



47r 



^ [{21 + 1) sin^ Si (A13) 

/=o 

-2{l + l)sin(5,sin(5,+iCos(J/+i - Si)] , 

which is a weighted integral over the differential scatter- 
ing cross-section. Estimating the maximum angular mo- 
mentum L that gives an important contribution to the 
sum, the authors of Ref. j37] approximated the above full 
expression for ax by assuming the phaseshift to be max- 
imal for I < L, and minimal for I > L, i.e. sin^ Si = 1 



and respectively. With this simplification, Eq. (A13) 
becomes 



47r 



(AM) 



which is in good agreement with the full numerical cal- 
culation [37]^ as long as £0 < e„ ^ 1; note that this 
condition is not satisfied in the Born regime (defined by 
large velocities, 1, or small couplings, :s> 1/2). 

For the models we consider in this work, we are par- 
ticularly interested in very small velocities (ct, <C e^). In 
this case the cross-section actually becomes velocity in- 
dependent and the above approximation fails (unless we 
are near a resonance, see below): when kR ^ 1 is satis- 
fied {R = ra^^ being the effective range of the Yukawa 
potential), the phaseshift is given by 



Si (y: k 



21 + 1 



(A15) 



^ In [191 138| the DM self-scattering rate was calculated in a dif- 
ferent way, using known expressions for or to estimate the ion 



drag force in plasmas |65l I66| . The results of |191 138| agree with 
|37l in the velocity regime that was of interest in these papers. 



All phaseshifts with I ^ are thus negligible compared 
to Sq and we have ctt — ci^o- The velocity independent 
effective potential cross-section is given by |67| 



crP°* = ai^o{kR<t: 1) = ina^ 



(A16) 



where we used sin^ Si ~ 5f = a 
called the scattering length. 

As in the case of WIMP annihilation, resonances can 
occur for low velocity scattering. Since the Schrodinger 
equation for / — scattering is the same as for annihila- 
tion, the resonances will occur at the same values of e^. 
Assuming these quasi-bound states appear at an energy 
Eq and have a width F, we are in the neighborhood of a 
resonance for \E — Eq\ ^ F, with E = m^v^ the kinetic 
energy of the system. In this region, the approximation 
(5o — ak as given by Eq. (A15) and Eq. (A16) receives a 



considerable contribution from the resonance: 



= ^0 + arctan 



F 

E{) — E 



The full transfer cross-section is thus given by 



pot 



where the resonance scattering cross-section [6 7) 

rcs^ 4^ r2-2afcF(J^-J?o) 
k^ {E- EqY + F2 



(A17) 



(A18) 



(A19) 



is no longer negligible compared to (tP°*, and should be 
taken into account. In fact, Eq. (A19) scales like l/fc^ 
and we can safely neglect ctP°* near a resonance. Exactly 
on resonance, we have 



a^'^^E = Eo) = 



47r 



(A20) 



which is the same as Eq. (A14) with L = 



Far away from a resonance we cannot use Eq. (A14) 



for ctt in the full velocity regime, since crP°* becomes 
velocity independent as u 0. In our analysis, we there- 



fore use Eq. (A14| down to the velocity at which err 



reaches the value of the asymptotic numerical solution 
for a-i=Q{v — )■ 0), beyond which we keep ctt constant as in 
Eq. (A16). Note that around the transition point, at in- 



termediate velocities, the full numerical solution actually 
results in ar being significantly larger than ai^o{v — 0) 
- both because of the behavior of (T;=o and because of 
contributions from ^ > 0; our approach is thus rather 
conservative in that it uses a lower limit on ar- 

For completeness, we note that a can actually also ex- 
hibit 'anti-resonances' for some values of v, and even 
(almost) completely vanish for very small velocities. 
This so-called Ramsauer-Townsend effect [33] can be ex- 
plained by a destructive interference between the poten- 
tial and resonance scattering amplitude, c.f. Eq. (A18l; 
indeed, it follows from Eq. (A19) that the resonance 
contribution is neg ative for T/JE - Eq) < 2ak. The 
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cross-section even disappears when the full phaseshift in 
Eq. (A17| becomes zero for T/{E — Eq) = tan(aA;). Since 



this only happens for special combinations of the param- 
eters that appear in the Schrodinger equation, we do not 
take this effect into account in this work. 



Appendix B: Constraints on new light particles 
coupling to leptons 

In this Appendix we summarize constraints on new 
light bosons i^, in particular those of the type that appear 
in our toy model introduced in Section flV A[ 

The strongest indirect constraints come from boson 
loops contributing to the anomalous magnetic moment 
of the muon, = (g^ — 2)/2. The deviation between 
the experimentally measured and theoretically expected 
value is Aa^ = a^^^P - af^ = (255 ± 63 ± 49) x 10"", 
which corresponds to a discrepancy of 3.2 times the esti- 
mated la error [4Qj . We will require that the contribution 
of a new theory does not worsen this discrepancy beyond 
the 5cr level, i.e. -1.45 x IQ-^^ < <™ < 6.55 x 10"^. 

The contribution of the pseudoscalar field to the mag- 
netic moment of the muon is given by |68] 



p ^ 1/2 Jk^ 

^/i cm 



dx- 



47r3/2 7q x2 + (1 - x)mj,/; 
and the contribution of the scalar field is given by 



1/2 Jk 



s2 



47r3/2 



dx 



2x^ 



-|- (1 — x)fn1/inn? 



(Bl) 



(B2) 



For a scalar of mass 0.1 (1, 10) GeV, the resulting limit 
on the coupling to muons is gf < 1.8 x 10^^ (6.7 x 10"'^, 
4.4 X 10~^); for pseudo-scalars, the limit is stronger by a 
factor of roughly 2 (for similar constraints, albeit usu- 
ally often considered for masses below 100 MeV, see 
e.g. Refs. [SSHTT] and references therein). Let us stress 
that these constraints only apply in the limit where the 
respective other particle's contribution to can be ne- 
glected. Since the two contributions have opposite signs, 
however, in principle even considerably larger values of 
gf could be allowed for suitable choices of rup and 5^ if 
one is willing to accept the necessary fine-tuning. In the 
case of exactly degenerate masses and couplings, e.g., the 
limits weaken to gf = g^ 0.003 (0.016, 0.15) for masses 
of = 0.1 (1, 10) GeV. 

The most important direct constraints on new light 
bosons derive from beam dump experiments, where the 
incoming electrons could radiate such particles and one 
tries to spot their decay products behind the stopped 
electron beam (see, e.g., Ref. ^ for an overview and 
Refs. [5ni[73] for a recent discussion). The only such ex- 
periment that currently can probe scalar particles heav- 
ier than around 100 MeV, however, is E137 at Fermilab 
[7Ij : for 100 MeV < < 400 MeV, couplings to elec- 
trons are excluded in a rather narrow band of roughly 



g'^ ^ 10~^ — 10~^. Future beam dump experiments may 
extend such limits up to masses of a few GeV and close 
the gap to the limits obtained by the muon g — 2 con- 
straints discussed above [73] [75J [75]. New constraints 
in this region may possibly also be obtained with a low 
energy electron-proton collider [69) . 

For masses > 2m^ ~ 210 MeV, the currently 
strongest constraints (apart from the small E137 win- 
dow mentioned above) derive from the BaBar search for 
T decays into light (pseudo) scalar particles that subse- 
quently decay into muons, e'^e" — >■ T(3S') — > 7^ and 
(f) — > fJ,^fJ.~ [77]. Since the final states are identical, 
the BaBar limits can be interpreted as limits on the 
direct production of light bosons in e'^e" ^cf); as- 
suming identical couplings to electrons and muons, this 
leads to roughly gf < 10"'^ - with, however, considerably 
weaker constraints for masses around the p resonance at 
« 770 MeV [75]. 

While our toy model is only loosely motivated by the 
excess |31| in cosmic ray leptons, let us mention that 
additional constraints in principle arise if one takes this 
connection serious and requires that the positrons (and 
electrons) from DM annihilation do fit the cosmic ray 
data. The fact that the cosmic ray antiproton spectrum 
[75] is consistent with the expectation for the astrophys- 
ical background [79[ puts severe limits on hadronic de- 
cay modes of the new light particles [HJ [SO] compared to 
those required for the leptonic modes (see, e.g., Ref. [42]): 
as a consequence, the coupling to quarks should very 
roughly be suppressed as g^ < 0.3 gf. The large annihi- 
lation rates that are required to fit the cosmic ray lepton 
data are also potentially in conflict with both gamma-ray 
and radio observations towards the galactic center [81) . 
gamma rays in the galactic halo due to inverse comp- 
ton scattering of the high-energy leptons [82j, as well as 
with measurements of the cosmic microwave background 
18, 62, 83, 84 or even big bang nucleosynthesis [85] . 

Completely independent of the cosmic ray connection, 
finally, very stringent constraints on low values of gi arise 
in principle from DM annihilation-induced gamma rays 
from galaxy clusters since in this case the smallest proto- 
halos form with very small masses and the annihilation 
flux from substructures is maximized; the corresponding 
constraints presented in Ref. [33], however, are rather 
model-dependent and not easily translated to our case - 
see the dedicated discussion in Section IVj 



Appendix C: Matrix elements and cross-sections 



In this Appendix, we present relevant interaction rates 



for the leptophilic model introduced in Section IV A 



The full expressions for the DM annihilation cross- 
section, which we used in our calculations, are somewhat 
lengthy and not very illuminating. Here, we therefore 
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only state the result up to O ^w^, rrv^/r 
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The fact that the two last cross-sections vanish for small 
velocities is simply a reflection of parity conservation. 
Note also that the direct s-channel annihilation into SM 
particles is strongly suppressed with {gl^^Y ^^^d thus neg- 
ligible. 

The matrix element for the elastic scattering of DM 
with standard model particles, if mediated only by a 
scalar, is given by 



\M'' 



^'^ (4m2 - t) {Amj - t) . 



(t — mIY 

In the case of a pseudo-scalar mediator, we find 



\MP\' = 



{t - mlf 



t^ 



(C4) 



(C5) 



In order to calculate the collision term for scattering pro- 
cesses to a sufficient accuracy, one only needs to evaluate 
the matrix elements at t = 0, c.f. Eq. (21 1. Therefore, 



only the scalar mediator gives a non-negligible contribu- 
tion to the scattering processes we are interested in here. 

The elastic scattering of DM with the light mediator 
particles (p can actually be much more efficient than the 
scattering with standard model particles - at least as long 
as rirj, is not yet strongly Boltzmann suppressed. This can 
directly be seen from a comparison of Eq. ( C4 ) with the 
matrix element for x4>s ^ x4>s which, in the relevant 
limit, is given by 



\M\U 



128 g 



4 



(C6) 



where uj is the energy of (ps- For the scattering with pseu- 
doscalar particles, x4'p ^ X^p; this expression becomes 



\M\ 



t=o 



128 gP^TO^w* 
nip — 4to^cj2) 



(C7) 



The inelastic scattering x4>p ^ x4>s is generally sup- 
pressed with respect to the dominant of the above two 



scattering modes; for degenerate masses 
e.g., we find 



\M\U 

s—rn'^--\-2'm^uj-\-rn\ 



32 .gfgfmX (i,^ - ml) 



Am^LO' 



(C8) 

In order to apply our formalism even to scattering with 
the light mediator particles, we finally need to make sure 
that these are still in thermal equilibrium with the heat 
bath. Unlike for the (assumedly) stable DM particles, 
equilibrium is in this case most efficiently maintained by 
(inverse) decay processes (j) -H- ££ and the first moment of 
the Boltzmann equation reads^ 



3iJn^ 



(27r)32w J (27r)32w J {2T:f2E 
x{2nf5(^\p-k-k)\M\l^u 

{[1-.9+H] [l-5+(^)] /(P) 

-5+H.9+(^)[l + /(p)]} (C9) 

-(r)(n^-n;'^) . (CIO) 



Up to one remaining integral in a;+ = wi -I- 0^2, the above 
phase-space integrals can be fully performed analytically 
and the thermally averaged decay rate is given by 



(Pk(Pk 5{E-uj-u) 



8(27r) n"" J J ujCo{uj + Cj) (eT + l) + I 



r^,Tmj 



duj4 



log 



cosh 



(Cll) 



4T 



cosh ■ 



4T 



(C12) 



In the above expressions, E = y (k + k)^ + m^, uj™ = 

^1 — Amf/ml^Lo]^ — and is the equilibrium 

number density which in the non-relativistic limit be- 
comes n'^ = m^K2{x)/{27T'^x). For very small tempera- 
tures, T ^ 7Ti0, (F) of course simply reduces to the decay 
rate at rest, for (j) ^ ££, given by 

3 



87r 

p2 



1 



(C13) 
(C14) 



We also calculated the annihilation rates and found that these 
are, as expected, much smaller than the decay rates. In fact, 
the tf> particles would in some cases never have been in thermal 
equilibrium if only 2 •<->■ 2 processes were taken into account. 
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